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An approach for particle-hole correlation functions, based on the so-called SCRPA, is developed. 
This leads to a fully self-consistent RPA-like theory which satisfies the /-sum rule and several 
other theorems. As a first step, a simpler self-consistent approach, the renormalized RPA, is 
solved numerically in the one-dimensional Hubbard model. The charge and the longitudinal spin 
susceptibility, the momentum distribution and several ground state properties are calculated and 
compared with the exact results. Especially at half filling, our approach provides quite promising 
results and matches the exact behaviour apart from a general prefactor. The strong coupling limit 
of our approach can be described analytically. 
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I. INTRODUCTION 

The advent of high T c super-conductivity, which remains unexplained theoretically in its essence, has spurred 
an enormous quest for developing many-body approaches capable to describe strongly correlated Fermion systems. 
Various formalisms have been applied in the past, each which its strengths and deficiencies (for a review see refJH). 

However, contrary to standard mean-field theory which is a commonly accepted lowest-order many-body approach, 
for correlation functions, such a generic method is missing so far. In this respect, any new and promising vistas are 
worthwhile to be pursued and elaborated. Indeed, in the recent past, a theory for two-body correlation functions has 
been developed bearing the characteristics of a generalization of Hartree-Fock theory to two-bodw dusters. In its roots 
this theory goes back rather far in time and has been promoted independently by several groupsBij. In the literature, 
it is known under various aaaies such as Self-Consistent RPA (SCRPA), Cluster Mean Field (CMF) and Equation 
of Motion Method (EOM)B~0. In itself it is an approximation to the so-called Dyson Equation Approach(DEA) to 
correlation functions where one replaces the full mass operator by its instantaneous part. However, only recently this 
method has found the attention and formal developments it deserves with, indeed, very promising results. The most 
outstanding of those is certainly the exact reproduction of the lowest spin-wave excitation spectrum, ujk = ^\smk\, 
known from the Bethe ansatz, of the antiferromagnetic Heisenberg chainEl. Moreover, also some simpler models 
have been treated successfully in this approachQ. Encouraged by these results, we thought it worthwhile to apply 
the method to the strongly correlated electron problem within the single-band Hubbard Hamiltonian with on-site 
repulsion U. 

Since the approach, which we hitherto want to call SCRPA, is based on non-linear equations for non-local two- 
body correlation functions, one understands that it is numerically very demanding. We therefore choose as a first 
application the one-dimensional Hubbard model for several reasons: 

• The exact solution of the ground state is known from the Bethe ansatz!. Therefore, a direct comparison of the 
SCRPA results is possible. 

• The numerical effort in one dimension may be expected to be more modest than in higher dimensions. 

• The experience gained in the Id case may help us to attack higher dimensions in the future. 
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The price to pay for this strategy is that one-dimensional problems are notoriously difficult to describe because of 
their extreme quantum character. As our method is not specifically designed for one dimension, we cannot expect it 
to reproduce particularities, such as Luttinger liquid behaviour. 

As we will see, our approach nonetheless permits to obtain interesting results in one dimension. They should, 
however, be judged in the light that in this first application to the Hubbard model we did not solve the SCRPA 
equations in full but applied a rather obvious and from the numerical point of view very simplifying approximation. 
Nevertheless, this approximation, known in the literature under the name of renormalizcd RPA, keeps the essentials 
of the self-consistency aspects. 

We demonstrate in this paper that the formalism allows for the self-consistent solution of a fully closed system of 
non-linear equations for two-body correlation functions. Moreover, important formal theorems are respected. Among 
those, we, for instance, cite the fulfillment of the /-sum rule (energy weighted sum rule) and of the Luttinger theorem. 

Other interesting results concern the strong coupling regime of the half filled chain. For example, the self- 
consistently calculated momentum distribution can be found analytically in the large U limit. It obeys n& oc cos/c 
with a proportionality factor of 4/J7 instead of 8 In 2/U of the exact result, known from the large U expansion of the 
Bethe ansatz solution, resulting in an error smaller than 30%. This result is the more astonishing as it was obtained 
with the renormalized RPA approach. It can be expected that it will be substantially improved once the full SCRPA 
solution is available. 

The paper is organized as follows. In section O, in order to make our paper self-contained, we will give a brief 
overview of the SCRPA theory starting from the equation of motion for a completely general Green's function. In 
particular, the explicit form of the self-consistent and renormalized R PA particle-hole propagators is derived in terms 
of a closed system of non-linear self-consistent equations. In section III, our approach is applied to charge and spin 
correlation functions in the Hubbard model. In section IV, we solve numerically the set of self-consistent equations 
for the renormalized RPA response functions for different fillings and interaction strengths. Our results are compared 
with the Bethe ansatz solution and with Quantum Monte Carlo calculations. For half filling and large U, analytic 
expressions are given for the momentum distribution function and the susceptibilities of our theory. In section [v|, 
we draw some conclusions and give an outlook on some improvements which are planned to be implemented in our 
approach. In appendix |A|, we outline the connection between the SCRPA and a variational ansatz which minimizes 
the energy weighted sum rule. Appendix [b| provides the explicit expressions for the free particle-hole susceptibility in 
one dimension. In appendix |^, we show how the analytic expressions which we derived for the strong coupling limit 
of our theory at half filling, are a rigorous solution of the renormalized RPA equations. 



II. THE DYSON EQUATION APPROACH TO MANY-BODY GREEN'S FUNCTIONS 

In this section, we briefly want to review the Dyson Equation Approach (DEA) to correlation functionsi. The DEA 
is increasingly used in-,tjbc many-body community and has recently produced interesting results in various domains 
of many-body physicsoQ. 

Let us start with the definition of a general causal (time ordered) or retarded many-body Green's function at 
zero temperature and at equilibrium (the generalization to finite temperature, using the Matsubara technique, is 
straightforward) , 

G\ B {tX) = ({A(t); B{t'))) c 

:=-i(0\T e A(t)B(t')\0) (1) 
G^(M') = ({A(t); B(t'))) let 

:=-i0(t-O (0| [A®, B(t')]_ e \0) , (2) 

where |0) is the exact ground-state and T e Wick's time-ordering operator. 

Here, A(t) and B(t') are arbitrary operators built out of any number of annihilation and/or creation operators of 
Bosons or Fermions or mixtures of both. Usually A and B will depend on one or several indices, and the notation 
(( A; B )) has to be considered as a shorthand for the matrix Green's function (( A a ; Bp )) where a and j3 run over the 
whole set of quantum numbers. The operators A and B can also be spin operators or even more general operators 
such as multi-component operators A = (Ai, A2, ■ ■ .) where the single components are again chosen according to the 
problem in question. For the derivation of a Dyson equation, however, we will choose B = A + . The case B =/= A + 
needs further considerations which may be important for the derivation of integral equations for vertex functions. 

The time dependence of the operators is given in the Heisenberg picture, X(t) = e zHt X e~ tHt , where the Hamilto- 
nian H is also completely general. It may describe relativistic or non-relativistic Fermi, Bose or spin systems or any 
system for which a Hamilton operator exists. 
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At equilibrium the two time Green's functions ([!]) depend only on the time difference such that their Fourier 
transforms are only functions of one frequency. These are the quantities for which we want to derive a Dyson 
equation. As the derivation is the same for either causal or retarded Green's functions we will from now on omit the 
upper indices. 

Let us establish the equation of motion for G AB {t, t')- 

i^ t {{A(t);B(t')})=6(t-t')([A,B}_ e ) + {{[A,H](t);B(t')}) , (3) 
or, in energy representation, 

u{{A;B)) u = ([A,B]_ e ) + {{[A, H]; B)) u . (4) 
On the rhs of (||) the commutator 

Af:=([A,B]_ e ) (5) 

plays the role of a norm matrix, and (( [A, H]; B))^ is a general Green's function containing the interaction once 
explicitly. For fermion like operators A and B we will use the anti-commutator Green's function (e = — ) and for 
boson like (e.g. a product of an even number of fermion operators) the commutator Green's function (e = +). 
Under the assumption that the inverse of ((A; B))^ exists, an effective "Hamiltonian" TCab(u) can be defined as 

■H AB (u 1 ) = (([A,H];B)) u ((A; B))^ . (6) 

The equation of motion (^) can thus be transformed into a Dyson equation, 

lo {{A- B)) u = N + H AB {u>) ((A; B)) u . (7) 

We stress again that the products on the rhs. of eqs. (||) and (0) are understood to be matrix multiplications. 
As we do not yet know how to determine the effective Hamiltonian Hab(^), the solution of the Dyson eq. (0), 

{{A;B))=L - Hab{u)X AT, (8) 



remains for the moment completely formal. In order to derive a more explicit and useful expression for TIab{<-o), we 
insert the inverse of the formal solution (0) in (0): 



Hab(u) = <( [A, H]; B)) u AT 1 ho - Hab{u) 

= H\ B {u) H 1 a] b {u) (9) 

The first part, Ti} AB (uj) = (( [A, H]; B )) Af^^-u, can be obtained from the equation of motion for the higher Green's 
function (( [A, H}; B)) u which is set up from "the right" this time (i.e. deriving with respect to t' instead oft): 

(( [A, H]; B)) u u = ( [[A, H], B]_ f ) + « [A, H]\ [H, B] )) w (10) 

If we adopt B = A + , the second^art of the effective Hamiltonian, TL 1 ab {lo) = (( [A, H]; B)) uj AT~ 1 Hab(u), contains 
only n-line reducible contributionsE3, with n being the number of fermion operators in A. Further, it can be showntLEI 
that the sole function of T-L l \ B (ijj) is to cancel all reducible contributions of 1~1 1 ab (lo). As the double commutator 
( [[A, H], B]_ e ) has no reducible contributions, we just have to put an index "irreducible" on the Green's function 
on the rhs. of eq. (EOI) to obtain as the final expression for the effective Hamiltonian: 



n AB (u) = {< P, H], B]_ f ) + (( [A, H\; [H, B] ))% J 



n A c B + nTM (ii) 



We see that the effective Hamiltonian ([11]) splits up in a natural way in an instantaneous part and in a truly 
dynamic (resonant) part. The latter contains scattering processes leading to imaginary potentials and corresponding 
real ones with a frequency dependence. 



3 



To obtain a better understanding of the various terms contributing to the effective Hamiltonian, let us analyze 
eq. ( |IT|) for the well-known case of the single-particle propagator, that is A = a\ and B = a^,. Since later we want 
to restrict ourselves to a non-relativistic fermion system let us consider a typical Hamiltonian 

H = ^i 12 afa 2 + ^ ^ ^1234 afa^ a 4 a 3 (12) 

12 1234 

where a, a + are fermion destruction and creation operators. The matrix elements ti2 and V1234 = W1234 — ^1243 of the 
kinetic energy and the two-particle interaction, respectively, are expressed in an arbitrary single-particle basis which 
comprises for example quantum numbers for momentum, spin, isospin, and so on. 

The norm matrix (||) is thus given by Afw = §iy ■ In this case, the effective Hamiltonian is the sum of the 
single-particle energy and the full self-energy. The static part of the effective Hamiltonian, expressed by the double 

commutator / [[ai, H], af,] + ^, yields the Hartree-Fock or mean-field Hamiltonian. We thus have recovered an 

important piece of the single-particle Dyson equation. Working out the second part of the effective Hamiltonian in 
eq. ( p"T| ) yields the following 2p — Ih Green's function: 

j ^ W1234 (((a^ a 4 a 3 ) t ; (aj, a+, a 2 >) t > )Y" «4'3'2'i' ■ (13) 

234 
2' 3' 4' 

As mentioned above, in eq. (|l3|), all reducible contributions to the effective Hamiltonian are removed and we obtain 
the usual irreducible self-energy T,u>(uj) of the single-particle Dyson equation by putting an index "irreducible" on 
the 2p — Ih Green's function in eq. ([l0|). 

The same scenario remains valid if we take for A and B more complicated operators like e.g. the density operator 
ay ■ Again only the (lp — lh) irreducible parts of the 2p — 2h Green's function in eq. ( ^p|) contribute to the effective 
Hamiltonian. 



A. Self— Consistent Random Phase Approximation 



As discussed above, the effective Hamiltonian splits, also in the general case, in an instantaneous and an energy 
dependent part. The instantaneous part can be considered as a generalized Hartree-Fock Hamiltonian (see below). 
Therefore, as a first approximation, one can try to solve this "HF problem" , neglecting the resonant part of the 
effective Hamiltonian. As we will see later, this allows us to solve e.g. the two-body problem on the level of a 
Schrodinger-like equation for a single-frequency Green's function, in contrast to the Bethe-Salpeter case where a 
three-frequency Green's function has to be determined. This means that we can introduce two-particle states with 
shifted energies. Therefore, the consideration of the instantaneous part of the effective Hamiltonian can be understood 
as a direct generalization of the common single-particle HF approximation to the more-body or cluster case. In .the 
past, this has been called Cluster Mean Field (CMF)El or Self-Consistent Random Phase Approximation (SCRPApo. 
In the remainder of this paper, we will adopt SCRPA as shorthand for our approach, which, for the two-body case, 
can be connected to a variational principle (see appendix [A|) . 

In analogy to the single-particle Green's function we thus can define a generalized n-body mean-field propagator 
by substituting the instantaneous part of the effective Hamiltonian back into the formal solution of the Dyson 

equation (|h: 



[A-B)) sc =\u;-n s A c B !> Af 



-1 



w-([[A,H],B]_ e )M-A N 



(14) 



Usually it is possible, as we will illustrate in an example below, to close the system of equations in the following 
sense: For a full set of operators A and B and for a two-particle interaction all expectation values in Hjfg can be 
determined self-consistently from the Green's function (M) via the spectral theorem. For retarded Green's functions, 
the spectral theorem at temperature (3 = l/iksT) readsEI: 

{AB) C = (AB) - (A) (B) 
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1 f^MMBK 
ir / 1 -ee~P u 



£^°-I|da;Im«4;S))f 
o 



~ J dcj Im((A;B»f 

— OO 



(15) 



indicates that we 
(A) (B). This 



The first of the eqs. ( jl5| ) is the well-known fluctuation-dissipation theorem. The superscript "c 
calculate a correlated expectation value, i.e. fluctuations of (AB) around their classical mean vi 
expectation value is known as cumulant or, in the sense of Feynman graphs, connected averages. As indicated in 
eq. Jl5|), the Fermi function in the spectral theorem for correlation functions (AB) C reduces to a step function as 
T — > 0. The commutator expectation values, such as ( [A, B] ), depend only implicitly on temperature, since the 
Green's function occurring in the spectral theorem ( |l5| ) is temperature dependent. 

From now on, we will restrict ourselves to the T = case, leaving the finite temperature consideration to forthcoming 
investigations. 



B. Energy Weighted Sum Rule 

The well-known energy weighted sum rule or /-sum rule connects the imaginary part of the exact Green's function 
((A; B)) u to the expectation value of the double commutator ( [[A, H], B]_ e ). Sometimes it is possible to choose 
operators A and B such that the double commutator ( [[A, H], B]_ e ) can be evaluated analytically. In this case, the 
sum rule may be used as a rigorous check for any approximative Green's function. 

Let us recall briefly the derivation of the sum rule. We can compute ( [[A, H], B]_ e ) using the spectral theorem 
( |l5| ) for the higher Green's function (( [A, H]; B)) w , 

OO 

( [[A, H], B]_ e > = -I J do; Im (( [A, H]; B))? . (16) 

— OO 

Inserting the equation of .motion (Q) on the rhs. and supposing the norm matrix to be real, we find the well-known 
energy weighted sum rule&S: 

OO 

([[A, H], B]_ e ) = ~J foulm((AiB)C (17) 

— OO 

We will now show that the sum rule (111) also holds for the SCRPA Green's function (( A; B ))^ c . From eq. (O) we 

sc 

see that ((A; B})^ satisfies the equation of motion 

u{{A;B))l c =M + ([[A, H], B}_ e ) U' 1 ({A; B)ff (18) 

rather tha n M ). Again, the norm matrix and the double commutator on the rhs. are real and u> independent. Inserting 
( ^8| ) into (|17|) and applying the spectral theorem ( |l5| ) yields the norm matrix on the rhs. (which is cancelled by A/" -1 ). 

From the above we see that, because of the double commutator structure of the effective Hamiltonian Tt sc , the 
SCRPA Green's function fulfills the energy weighted sum rule ( |l7| ) practically by construction. 



C. Particle— hole propagator 

As a concrete example, we will derive the SCRPA expression for a particle-hole Green's function a+ <Xk) a-k> a p' )Y l ^ t 
in a fermionic system with general two-body interactions as described by the Hamiltonian ( |i"2| ) . Since the operators 
defining the Green's function are pairs of fermions, we will use the commutator Green's function (e = +1). 
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For a homogeneous system, the kinetic energy is diagonal in momentum space, 

tkk> — 5kk' £k , 

with k standing for momentum and other quantum numbers such as spin. 
The norm matrix is also diagonal, 

Nkpk<p> = ^ [dp a,k, a£, Op/] _ 
= Skk/ 5 PP > (n p - n k ) 



(19) 



where 



>ik 



< a t «fc ) 



(20) 



(21) 



stands for the occupation numbers. The effective SCRPA Hamiltonian, introduced in eq. (|llj), can be worked out for 
the Hamiltonian (|12|). Using summation convention, this yields 



H 



kpk'p' = ( H], a+ ctp/]_ ^ (ry - rife/) 1 

= <5fcfc' <5 PP ' (efe — e P ) + (rip — rife) v P 'kk' P 

— Spp' Vkq 2 q 3 q4 { a k> a q~ 2 a l3 a <j4 ) + 2 ^91 92 93 P ( a ^ a ^ a 93 a p' ) 



1_ / _|_ 1 \c 1_ / 1 

2 Vfe?4 \ a k' a p a 93 a ?4 / + 2 U 9l92j>fc' \ " 



a qi a q ^a k a p > ) 



v kqiq 3 k' { &p a q~ 2 O93 a p' ) v p' qnqzP { a k' fl 92 °93 a fc ) 

where denote the Hartree-Fock corrected single-particle energies, 



rife/ 



£fe + Wfcg/jg Tig 



(22) 



(23) 



As we will see in section [ID , the second term in eq. fl22[), (n p — rife) Vp'fcfc'p, will lead us to a RPA-like theory. 
We will use the term RPA in a slightly broader sense than usual and already account for the exchange term of the 
interaction, see ecu- (^3|). The term in brackets, in contrast, contains exclusively correlated expectation values (i.e. 
cumulant averageala), 

( a ti a t a q3 a 94 ) C = ( a+a+a q3 a qi ) - [( a+a qi ) ( a+a q3 ) - ( a+a q3 ) ( a+a qi )] 

not taken into account by usual RPA-like theories. 
The SCRPA Green's function 



(24) 



Skpk'p'i^) = (( a t a k\ a t> a P ' )) 



sc 



(25) 



u ~ n kpk'p' 



Once it is determined, all elements of 



defined in eq. (^) , can now be obtained by inverting the matrix 

the effective Hamiltonian (g2 ) and the norm matrix (|2(]) can be calculated via the spectral theorem ([L5|). Moreover, 
it is possible to derive an explicit expression for the occupation numbers n p by summing the diagonal elements of the 
norm matrix, n p — rife, over the index k. 

With the commutator spectral theorem (|l5|) and the fc-space volume 



k 

we get for the occupation numbers 

00 



(26) 



(27) 



G 



In a continuous system, it is necessary to introduce a cutoff in order to keep the fc-space volume V finite. In lattice 



systems, as will be seen in section III, V is finite, since the fc-sum is restricted to the first Brillouin zone. 

The correlation functions ( a^ 1 a^ 2 a q3 a qi )° in eq. ( |22] ) are connected to those, which are accessible via the spectral 
theorem (151), 



( ( a tl °94 ) • (at a 93 ) ) C = ( a tl a i4 a t ) - ( a t a Q4 ) ( °?3 ) • (28) 

This yields 

oo 

( a i a i a 53 a 94 ) C = ~ J dw Im (( a+ a Q4 ; a+ a g3 - ( a+a Q3 ) [J g2(Z4 - ( a+a 94 )] , (29) 
o 

where for a homogeneous system the last term on the rhs. is related to the occupation numbers by momentum 
conservation. 

The system of equations is now closed and can be iterated to self-consistency. We therefore start with an assumption 
for the expectation values in H sc and A/". The SCRPA Green's function, obtained by matrix inversion, then allows 
us to calculate new values for the elements of the effective Hamiltonian and the norm matrix by applying the spectral 
eqs. ( prb and (p9|) . The new Tt sc and Af lead us to the next approximation for the Green's function and so forth. 



D. RPA and Self-Consistent RPA 

We now analyze the contributions to the SCRPA Green's function by rewriting eq. (0) as an integral equation: 



k 1P1 
»=2P2 



where Qkpk'p'i^) nas the structure of a free particle-hole Green's function, 



Gkpk'p'i^) = S PP > j-^- " fc , . n+ ■ (31) 

y p u) - (efc - e p ) + «0+ 



The integral kernel I^t?p 1 k 2 p2 re P resen ts the interaction occurring in the effective Hamiltonian (|2S 

K-kpk'p' = i n P - n k) 1 ['Hkpk'p' ~ hk' &pp' (efc - e P ) ] ■ (32) 
Since the kernel JC splits up into a RPA-like part 

fclepk'p' = Vp'kk'p (33) 

and a remainder, /C c , which only contains correlated two-body densities, it is convenient to rewrite eq. (^) in a 
different way (using matrix notation) 

gSC = qKPA + qKPA ^SC (34) 
gRPA = gO + g0]C KPAgKPA ( 35 ) 

At this point, we emphasize that the eq. (|3^) has exactly the same structure as the usual RPA equation. However, in 
our theory the occupation numbers can and will, even at zero temperature, be different from the Hartree-Fock values 



,HF 



= 0(E F - e k ) . (36) 



In the following, we will consider the eq. (|35| ) as generic for the RPA whatever the occupation numbers will be. We 
will label it pure RPA if the occupation numbers are fixed to their HF values (|36|). In contrast, a theory in which the 
occupation numbers are determined self-cansistently from the RPA Green's function or contain correlations in any 
other way will be called renormalized RPAu. 

In SCRPA, eq . (pj ) is coupled to the RPA eq. (|^), upgrading the (renormalized) RPA to the Self-Consistent RPA. 
Comparing eq. (p0|) to eq. (35) clearly shows that the RPA structure of the solution is preserved when passing over 
to self-consistent RPA: 
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e sc H = [i - gV)£ s T s°n 

g RPA (oj) = [i - g°(uj) k rpa ] - 1 g°(uj) (37) 

The improvement contained in the integral kernel K, c can be interpreted by the following connected diagrams: 





(38) 



The hatched circle represents a correlated particle-hole propagator and the dot stands for the antisymmetrized 
interaction v. The crosses represent Kronecker symbols in momentum and other quantum numbers, and <5-functions 
in time. We see that the first graph in eq. (|3^) corresponds to a coupling of the single-particle motion to the density 
fluctuations, i.e. a self-energy correction, whereas the second graph describes an induced (screened) interaction. Of 
course analogous graphs exist where the interaction is attached to the hole line. Again, we want to emphasize that 
all terms in eq. ( pq ) are instantaneous. We obtain the second-order contributions in replacing the hatched circle in 
eq. (p8|) by an interaction dot: 




(39) 



Cutting the graphs in eq. (|39j) between the two interactions occurring one an infinitesimal time after the other 
illustrates the instantaneous coupling of the lplh and 2p2h spaces. Solving the eqs. ( |35| ) and ( |3~4| ) self-consistently 
thus constitutes a partial resummation of the interaction to a very high order. 

Just as the Hartree-Fock self-energy for a single particle can be constructed from a two-particle interaction term 
by attaching an outgoing to an incoming line, viz. 





(40) 



we may interpret ecu (|3§|) as the Hartree-Fock field for density-density fluctuations (this point of view has actually 
been adopted in ref.B) . In analogy, we can reconstruct the loops in eq. d38| ) by closing two density fluctuation lines in 
the following first-order terms for the interaction (which can be obtained from perturbation theory): 





(41) 



Considering all exchange terms it is possible to reconstruct exactly the effective Hamiltonian (|22|) which therefore 
represents the mean-field Hamiltonian of a gas of quantal fluctuations present in any many fermion system. 



III. APPLICATION TO THE HUBBARD MODEL 



In this section, we will apply the SCRPA, developed in section |n[ to density-density correlation functions in the 
Hubbard model. I— . 
The single-band Hubbard Hamiltonian describes electrons hopping on a lattice with an on-site interaction C/li-a: 
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H = XX' "«'"■"' + ^X^T^l ( 42 ) 

ija i 

where af a and denote the creation and destruction operators for an electron with spin a on site i, respectively. 
The occupation number operator on site i is defined as hi a = af a ai a - In the following, we will restrict ourselves to 
nearest neighbour hopping 

Uj = -t(Sj ti+ i +5j ji -i) , (43) 

repulsive interactions U > 0, and zero temperature. We will work with % = 1, set the lattice spacing to unity (a = 1) 
and measure energies in units of the hopping integral (t = 1). 
After Fourier transformation, the Hamiltonian reads 

H = ^2e k a+ a a ka + — ^ a£ T a k+qT a+ ap.qj . (44) 

ka" kpq 

Notice that N is the number of sites (ions, not electrons) and all momentum sums run over the first Brillouin zone 
unless indicated differently. The single-particle dispersion relation in the hyper-cubic lattice is given by the Fourier 
transform of the nearest neighbour hopping matrix element ([43]), 

d 

£k = — 2^cosfci , (45) 

i=l 

where d denotes the dimension. 



A. Charge— and longitudinal spin— density correlations 

In the following, we will examine the behaviour of charge- and spin-density fluctuations in the Hubbard model. 
Therefore, we will introduce the density operator 

p qCT = a k t r a k+q CT , (46) 

k 

which is the Fourier transform of the Wannier number operator h icr . It will be used to describe charge and longitudinal 
spin fluctuations. Summing over the spins gives rise to the charge susceptibility 

x ch (q, W ) = i (( ( PqT + p qi ); (p+ + p+) . (47) 

As from now on we will use only retarded Green's functions, we will omit the superscript "ret" on correlation functions. 

The z-component of the spin on site i can be expressed as the difference between the number of "["-spins and .[-spins 
on that site, or, after Fourier transformation, 

- | (Pq T " Pql) ■ (48) 

Correlations between the z-components of the spins are described by the longitudinal spin susceptibility, 

X sp (q,-) = ^«5-^ + L 

We will now examine the charge and longitudinal spin susceptibilities, x ch (q, uj) and x sp (ci; w ) m SCRPA. Therefore, 
we introduce the particle-hole Green's function 

£k<x P <x'(q,^) = (( a^flk+q^; ap+ q(J < a P a' )) • (50) 



In contrast to the general particle-hole propagator defined in section II C , we will now account for momentum 
conserv ation right from the beginning. As the derivation of the SCRPA propagator is completely analogous to 
section 



II C, we will only state the results. 
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The norm matrix (EG) is given by 



A/kff p<7 ' (q) = 8 kp 5 aa > (n k <r - n k +qcr) 



(51) 



where Uko- denotes the occupation number for the Bloch state k with spin a defined in analogy to eq. (^lj) . The effective 
Hamiltonian, which was defined for a general two-body interaction in eq. reads for the Hubbard Hamiltonian 



Kkapa'il) = <5kp<W [£k+q - £k] + <5<7,-tr' (rikff - ™k+qcr) 
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N 



-<W<W — ( (a 



N 



"k+q-q'o^ktr + a k+qo- a k+q'o' I Pq-q',-<7 



J7 



+ \ c 

I p+qer a k+q<T Pp-k,-<T f 



u 

N 



{ ( a ia a ^ 



q-q'cr a k+q'o'"k+qo' 



ak4 



p+q-q',-cr a P^- cr a p+q,-<r a P+q',-a~ 



(n po 



'p+qer' 



(52) 



with pqa- being the density operator introduced in eq. (^6|). The Hartree-Fock corrections to the single-particle 
energies cancel because of the on-site (and thus momentum independent) interaction. 
The spectral theorem yields for the occupation numbers (see eq. j27|)) 



oo 

~~ ^n^2 J dw Im ^k<Tk<T(q,^) , 



(53) 



where ( n a ) denotes the number of c-electrons per site. In the paramagnetic phase, spin-broken expectation values 
like ^ak| ak l^ vanish, and we obtain from equation (^9|) for the correlation functions occurring in the effective 
Hamiltonian (K 



a kT a k+qT a p+qi a Pi 



dui Im5ktpi(q,w) 



(54) 



As in section II C, the system of equations is now closed and can be iterated to self-consistency, since, on one hand, 
we are able to calculate all elements of the effective Hamiltonian W k ^ po ./(q) from the particle-hole propagator, and, 
on the other hand, the latter by a matrix inversion for every q and lu: 



£w PCT '(q^) = [ w -^ pff '(q)] 1 (v -n p+Ci a>) 



(55) 



We will see in section III C that the corresponding system of equations for spin-density correlations is not closed onto 
itself, but couples bac k to t he charge-density correlations. 

As shown in section [IC (see eqs. ( 3^j35| )), it is advised to first calculate the RPA particle-hole propagator before 
solving the full SCRPA problem ([55). For the Hubbard interaction, the RPA kernel /C RPA defined in eq. ( pjj ) is 
nothing but the interaction per site 



/cL p pv(q) = 5 a ,. a Z 



Eq. ( p5| ) can therefore be written as an integral equation coupling Gy^ pcr {<h w) and G^~^ p<7 (q> w ) : 



(56) 



(57) 



where Q^ a (q,oj) defines the renormalized free particle-hole propagator 



<?L(q,") = 



(n kfT - n k+qg ) 
w - [£k+ q - £k] + «0 H 



(58) 
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In terms of Feynman graphs, this means substituting the SCRPA kernel , (q) in the integral eq. (| 

expression, which is nothing but a spin-flip interaction, represented by a dot, 



by its RPA 



k, cr 



k,<r 



u 



P, a 




q, a 



p + q, a 



k + q, a 



q : cr 



p, a 



(59) 



The physical interpretation of eq. (|57|) is that in RPA the Hubbard interaction flips the spin of the electron 
on each and every scattering process. This is certainly a good approximation for an electron propagating in an 
antiferromagnetically ordered state. In this case, we see from Fig. [j] that a J, -electron added to the "["-electron on site 
i cannot hop off because all their neighbours have spin J,, too. Therefore it must be the "[-electron which hops to a 
neighbouring site. Arriving there, it can only hop the same way back to its original site. Otherwise, the "[" -electron is 
surrounded by other "["-electrons and thus frozen in. As this process continues, an extra electron propagating in an 
antiferromagnet from site i to site j flips every spin on its trajectory. Thus, neglecting higher-ordet loop trajectories, 
the electron's path is completely retraceable. This case was first examined by Brinkman and Ricaij who showed that 
this "retraceable path approximation" is accurate for walks up to length twelve for the analogous case of an extra hole 
propagating in an antifcrromagnetic spin configuration. Moreover, they showed that even if the spins are randomly 
distributed rather than antiferromagnetically ordered, the dominant contribution to the hole Green's function comes 
from retraceable paths. 

In this line of reasoning, the dimensionality of the system plasma a crucial role. In one dimension, antiferromagnetic 
long-range order is forbidden by the Mermin-Wagner theorernEj. Nevertheless, as there are no loop trajectories in 
one dimension, the retraceable path approximation becomes exact for any spin configuration. 

In dimensions d > 3, antiferromagnetic ordering is possible. As the number of nearest neighbors increases, loop 
trajectories become less probable. Therefore, the retraceable path approximation gets exact to order C(l/d 4 )tZl. 

Moreover, in higher dimensions the correlations are weaker than in lower ones. As we derived the RPA kernel 
K-ka pa' (q) by neglecting the correlations present in the SCRPA, we expect it to be more accurate in higher dimensions 



than in lower ones. We will see in section IV that, even in one dimension, the renormalized RPA solution shows a Mott- 



Hubbard transition at a finite critical interaction strength, which, for half filling is of the order of the bandwidth. 
This scenario, which can be considered as generic for the RPA, is certainly wrong for the one dimensional case. 
Nevertheless, in higher dimensions i t j j oy ld be quite realistic. There, this viewpoint is also supported by methods like 
e.g. the Hubbard-III approximationllao. 

Iterating the integral eq. we can decouple the equations for £?kJp ff (q, w) and G^ap-cri^ u>), 



<5k P Gbi<i,v) 



N 



(60) 



where we introduced the renormalized non-interacting susceptibility 

x °(q,c) = i^C(q^). 



(61) 



Eq. fcG) can be solved explicitly for the particle-hole Green's function, yielding 



G^M,u) = GlMM 



5k P + —G p „{q,Lu) - 



1 

l-tfX?(q,w)tfx?(q," 



and 



(62) 



Finally, we have to determine the occup ation numbers nko- self-consistently from the RPA Green's function (p2|). As 
will be explained in more detail in section IV A, this is a somewhat delicate procedure. Indeed, initiating the iteration 
cycle, at small U , with the Fermi step for one inevitably picks up some spin instabilities when summing over q 
in eq. (|5^). These instabilities correspond to poles in the RPA Green's function at purely imaginary frequencies. This 
entails that we may create some unphysical values when applying the spectral theorem to such a Green's function. For 
example, integrating ImQ^^_ a (q,u)) over the whole w-axis yields an unphysical finite value if the Green's function 
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has poles at imaginary frequencies. As this integral is connected via the spectral theorem (15) to a commutator, we 
know that it has to be zero. When iterating the system of non-linear equations, it turns out that this pathology is 
cured. The occupation numbers get rounded thus weakening the interaction in such a way that finally the imaginary 
spin poles disappear. At self-consistency, the integral of Im <?kJp_ CT (q, over all frequencies vanishes for every q, as 
it is expected from the spectral theorem. 

For completeness, we state that fixing the occupation to their Hartree-Fock values brings us from renormalized 
back to pure RPA: 



HF 



e(£ F - e k ) (63) 



In section IIIB, an energy weighted sum rule is shown to be fulfilled in both, renormalized and Self-Cons istent 



RPA. In pure RPA it is only fulfilled as long as all cigcnfrcquencics arc real. As will be pointed out in section IV F , 
this is closely connected to the problem of RPA instabilities. 

If we restrict ourselves to the paramagnetic phase, we have — n^i implying X|(Qi w ) — xjCli w )' We can thus 
recover the usual RPA structure for the susceptibilities x ch (q, w) and x sp (q,uj) by combining the two equations j6^ ) 
and summing over k and p: 



x sp 



kper 

_ 2x°(q^) 
l-tfX°(q,w) 

kper 

^ M (64, 



1 + C/ X »(q,«) 



B. Energy weighted sum rule 



The fact that the density operators p qcr commute with the interaction term of the Hubbard Hamiltonian 



gives 



us the possibility to evaluate 



analytically. In analogy to section II B, we thus can establish an 



energy weighted sum rule for the exact susceptibilities x ch (q, w) and x sp (q, We will now discuss the fulfillment of 
the sum rule for the SCRPA, the renormalized RPA, and the pure RPA Green's functions. 
In the first place, we calculate the double commutator 



k 



[Ck+q — £k] (n>ka — «k+qc) 



8 aa > 2 ^2 ( C0S( 7j 

i=l 



ki 



(65) 



where ki are the components of the vector k. The only contributions to the double commutator come from the kinetic 
term of the Hamiltonian (|44|), since, as mentioned above, p q(T commutes with the interaction. Remark that for a 
k 2 /(2m) dispersion law the rhs. of (pq ) yields the well-known result ( n a ) q 2 jra (see ref£j). 

As pointed out in section p B for the general case, the expectation value of the double commutator (|65l) is related 
to the imaginary part of the exact retarded susceptibility by an energy weighted integral (see eq. (|TT 



du> ui Im 



Pqcr'j Pqcr' 



(66) 



The SCRPA Green's function was shown to satisfy the sum rule ( pq ) as well, since the double commutator on the lhs. 
can be expressed by the SCRPA Hamiltonian (52) and the norm matrix (51): 



12 



[p qCT , H], p { 



qcr' 



E E ^p lCTl (q)A/p lCTl p^(q) 
kp picri 

E W kap CT '(q) Op*' - «P+q<T') 

k P 



(67) 



In the view of our formalism this may seem evident. However, theories which generalize the RPA approach do not 
necessarily fulfill the /-sum rule. Above all, we notice that in eq. ( |67|) all terms containing correlations from the effective 
Hamiltonian 7ij^ p(T , (q) cancel when summing over k and p. In renormalized RPA, on the other hand, we neglect 
these correlations right from the beginning (see eqs. (|5(3| , |57j )). Consequently, the renormalized RPA susceptibilities 
fulfill the energy weighted sum rule (pfl), too. Moreover, it is well-knownll3E3 that the pure RPA susceptibilities satisfy 
the sum rule (|66|) if the expectation value on the lhs. is evaluated with the Hartree-Fock ground state wave function. 
However, this statement only holds true as long as all RPA frequencies are real. 

Finally, we find the energy weighted sum rules for x ch (q, oj) and X sp (qi w ) by combining eq. ( |66| ) for the different spin 
configurations according to the definitions of the charge susceptibility ( |47| ) and longitudinal spin susceptibility ( ptiT 
respectively: 



2E (oosgj-l) it 



i=l 
d 



doj oj Imx (q, oj) 



\ E (cos 9i - 1) ) = -- f doj oj Im X sp (q, W ) 



(68) 



On the lhs.. we introduced the shorthand 



# = 



T? E £ {^} n {fe.}^ 



(69) 



standing for the contribution to the mean kinetic energy per site provided by the electron motion in i direction. It 

should be stressed that / t l \ depends on the occupation numbers nicer and thus implicitly on the Green's function. 

In contrast, for the usual k 2 /(2m) dispersion law, the double commutator in eq. (65) depends on the mean number 

of electrons per site (i.e. the filling) instead of an d is therefore a model independent quantity. In our case of 

a cosine dispersion law, however, the q-dependence of the sum rule provides a check which does not depend on the 
Green's function or any other assumption. 



C. Transverse spin— density correlations 

In this section, for completeness, we will shortly discuss the transverse spin response. The spin-flip operators Sf 1 
may be substituted in the usual way by combination of an annihilation operator of a a electron and a creation operator 
of a —a electron. After Fourier transformation to momentum space, we obtain 

= E a kT ak +li aild = E a k+qi a lT ' ( 70 ) 

k k 

Correlations between the spin-flip operators give rise to the transverse spin susceptibility 

X + -(q,u) = b((S+;S- )) u . (71) 
In order to examine the transverse spin susceptibility in SCRPA, we introduce the following correlation functions: 

->trans 



(1. = (( < a k +q ,- ff ; a; +(i _ a , a pa , )) (72) 



In complete analogy to the charge-density case (see section III A ) , we calculate the norm matrix by evaluating the 
commutator of the two operators defining the Green's function (see eq. (P0|)): 

ACpMq) = *k P <W - n k+q ,- CT ) (73) 
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The effective SCRPA Hamiltonian reads in analogy to eq. (E 



WkJp S CT '(q) — SkpSaa' [ek+q,-CT — Cktr] — ^crcr' ("ker — «k+q,-cr) 



[7 

iV 



<5kp^(T(T' at^Z ( a k+q,-CT a k+q',-(T Pq-q'er + a k+q-q'er a k<x Pq-q',-<T 



q' 



q' 



Saa ' ( a kcr a k+q-q'cr - a£ +q , ^Ctk+q-t, 



u 



tier -a' ( a^ap-a S ] 



k P 



"p+qer 



a k+q.-cr S'p_ k ) 



a p+q-q'cr a pcr a p+qi _ a .a p + q ' i _ (r 



(^pcr' ^p+q,— <x') j 



(74) 



with p qcr denoting the density operator introduced in eq. (^) and iS q being a spin-flip operator defined in analogy 
to S+ (see eq. ((7C 



^q - ^ a k<j a k+q -<r 



(75) 



The Hartree-Fock corrected single-particle energies defined in eq. (|23|) are given by 

fker = e k + U (n_ CT > . 



(76) 



In contrast to the charge-density case, not all elements of the effective Hamiltonian (|74|) can be determined self- 
consistently from the Green's functions ((72j). The calculation of terms like ^ a^ak+q-q'T a p+ q - q '| a pT ^ > f° r example, 
cannot be performed with the spectral theorem for spin-flip Green's functions, since they contain always the same 
number of f a nd J spins. These terms can, however, be determined from the charge-density Green's function intro- 
duced section [II A . By this means, the transverse spin susceptibility is coupled to the charge-density susceptibilities 
X° h and x sp - I n this work, however, we will not further investigate the transverse spin excitations. 



IV. RESULTS FOR THE CHARGE AND THE LONGITUDINAL SPIN RESPONSE IN THE HUBBARD 

CHAIN 



As a first application of our general formalism, we will calculate the charge and longitudinal spin correlation 
functions X ch (<7; w ) an d X sp (<iS w ) m the one-dimensional Hubbard model. This will also serve as a test of whether our 
formalism is well behaved in a num erical sense. 

It was explained in section HI A that the first step will consist in calculating the Green's function Q kapu i (g, w), 
introduced in eq. (|50|), on the level of the renormalized RPA. In this paper, we will not go beyond this approximation. 
Indeed, the numerical solution of the full SCRPA problem turns out to be quite enormous and certainly needs a 
major numerical effort which is intended to be invested in future work. Nevertheless, the main characteristics of the 
self-consistency cycle are already present on the level of the renormalized RPA. 

We thus will determine self-consistently the RPA Green's function (^2|) together with the occupation numbers n^a ■ 
Due to the self-consistency, the occupation numbers are "renormalized" which modifies the Green's function G^aili w ) 
and the susceptibility x°(?,w) occurring in the RPA propagator (p^). 

In contrast, in pure RPA, i.e. if we fix the occupation numbers to their Hartree-Fock values given by eq (|6J 
Gka-il^) an d X°('?! a; ) are identical with the free particle-hole propagator and susceptibility, respectively. 

After a brief overview of the numerical method we will discuss the results for the renormalized RPA in the infinite 
Hubbard chain in the paramagnetic phase, i.e. rik] — n^. 



A. Numerical Method 

In order to determine the renormalized RPA Green's function, we have to solve the RPA eqs. ( p8| , |6l| , p2[ ) consistently 
with the equation for the occupation numbers, (|5^). 

Therefore, we start with an assumption for the initial occupation numbers and set up the following iteration cycle: 

• Set up the renormalized free Green's function ( ^8| ) from the current set of occupation numbers. 
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• Calculate the renormalized free susceptibility ( |6l| ) by integrating G^ail) w ) over Note that for the computation 
of Imx°(g,u;) it is convenient to introduce a small but finite imaginary part "i0 + " in the denominator of 

• Set up the RPA Green's function ( |62|) from Q^ :cr (q, u>) and X°(?, ui). 

• Compute a new set of occupation numbers by performing the spectral integral in eq. (p3|). This task is also 
simplified by assuming a finite "z0 + " . 

• Repeat the iteration cycle from the beginning until self-consistency is achieved. 

Usually, we will start with a small interaction, e.g. U = 1, and a Fermi step for the occupation numbers. For this 
U, the RPA equations are then iterated to self-consistency. The result is used to initialize the occupation numbers 
of an iteration cycle with a slightly higher interaction. We thus increase the interaction in small steps, iterating each 
time to self-consistency, until the desired value for U is reached. This procedure has the advantage that the spin 
response, in contrast to pure RPA, remains stable throughout the whole calculation. 

The momentum integrations are performed by summing over a grid of uniformly distributed points in the first 
Brillouin zone. The number of points is typically 256. 

Energy integrals are computed using a grid of points obeying a Lorentzian distribution peaked at ui — 0. The 
number of energy points is typically 2048 and the half width half maximum of the distribution is about twice the 
bandwidth, i.e. 8. 

As mentioned above, it is convenient for computational purposes to introduce a small but finite imaginary part 
"i0 + " in the denominator of the renormalized free Green's function G kcr (q, w). Typically we use values of the order of 
magnitude of i/16 which corresponds to 1/64 the bandwidth. It can be shown that within our resolution the results 
are not affected by slight variations of this value. 



B. Occupation numbers 

Fig.^ shows the momentum distribution function 

n k = \ (n feT + n kl ) (77) 

of the half-filled Hubbard chain. The renormalized RPA results are compared with Quantum Monte-Carlo calcu- 
lations. We see that, for small U, the renormalized RPA momentum distribution has a discontinuity at the Fermi 
edge. This is typical for Fermi liquids, and thus indicates a metallic behaviour. For large U, the renormalized RPA 
momentum distribution is continuous, and our theory predicts an insulating ground state. Increasing U in small steps 
(1/4), we deduce that within our numerical momentum resolution the discontinuity at kp vanishes at U ~ 3. We thus 
find a Mott metal-insulator transition at an interaction strength which is slightly smaller than the band width (i.e. 
4). This is in good agreement with approximations designed for the Hubbard model in higher dimensions, like e.g. 
the Hubbard-III solutionaHni3 In one dimension, however, the exact solutionis predicts an insulating ground state for 
any finite U. Consequently, the exact resjiLts. for n k , known from Quantum Monte-Carlo (QMC)EH for finite U and 
from Bethe ansatz in the limit of large £/E3'E3, show a smooth behaviour over the whole k range for all interactions. 
This disagreement should be judged in the light that the effective Hamiltonian of the renormalized RPA neglects 



all correlation functions. As discussed in section [II A, this approximation is expected to be far better in higher 



dimensions. Nevertheless, in the strong coupling limit, our theory reproduces the cosine-behaviour for n&, known 



from the Bethe ansatz expansion, apart from the prefactor (see section IV G). 

Away from half filling, the renormalized RPA predicts a similar scenario: For small U, the momentum distribution 
function shows a discontinuity at kp, which now persists up to higher interaction strengths as in the half filled 
case. For the quarter filled chain, e.g., this jump lasts up to U ~ 4, as can be seen from Fig. [| The discontinuity 
occurs precisely at the same k = kp as the step in the momentum distribution of the free Fermi gas. Therefore, the 
renormalized RPA satisfies the Luttinger theorem in one dimensioned. 

Again, the Mott-Hubbard transition predicted by our theory is in disagreement with the exact behaviour away 
from half filling, which is known to be a Luttinger liquid for all interaction strengths. The latter is characterized by 
a power law singularity in the momentum distribution at the Fermi points, which is expected in the strong coupling 
limit from both, Bethe ansatz expansionsEEfO and QMC calculations^]. For finite interaction strengths, it was also 
detected by QMC studies of infinite Hubbard chainsEfi. 

■Calculations of finite chains indicate a discontinuity at kp decreasing only very slowly with increasing chain length 
AotlJ. Therefore, they have great difficulties to detect the Luttinger liquid behaviour. There are certain similarities 
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between finite chains and calculations using a finite number of points in the Brillouin zone as e.g. the numerical 
solution of the renormalized RPA equations. Nevertheless, our calculations indicate a finite slope of the momentum 
distribution on both sides of the Fermi points, and we do not find any signature of Luttinger liquid behaviour. 



C. Renormalized free susceptibility x° {<!,<*>) 
1. Imaginary part 

Before discussing the RPA response functions x ch (9,w) and x sp (?,^) let us analyze the renormalization effects in 
the free susceptibility X°(<li u> ). 

In the (q, w)-plane, the imaginary part of x°(q, u1 ) is restricted to the region, where particle-hole excitations exist, 
i.e. there must be a k for which uj = £k+ q — £k is satisfied, or, 

M < 1 4 sin || . (78) 

If the occupation numbers are step functions like in the free Fermi gas, a second boundary condition is provided by 
the fact, that in the particle-hole continuum described by eq. (|7q) there must be at least one fc-vector for which the 
numerator of the free Green's function j5q), (jiha — ?ife+ga-)j does not vanish. This yields 

M > 2 |cos(fc F ) - cos(fc f - \q\)\ . (79) 



We have seen in section [II A that the Hartree-Fock corrections to the single-particle energies cancel. Hence, the 
denominator of the free Green's function (^8|) does not change throughout the renormalization process, and eq. J7§| ) 
represents a rigid boundary for the imaginary part of the free susceptibility. On the other hand, any rounding of the 
occupation numbers will directly affect the second boundary condition, (179) . This behaviour of the imaginary part 
is illustrated in Fig. ^ for half filling and q = n/2. We see that Imx free (<?, ui), represented by the dotted line, is only 
nonzero in t he dom ain in between the two boundaries given by eq. ([78]) and eq. ([79]). 

In section |lV A , we explained that for technical reasons we have to set "i0 + " in the denominator of Q® a (q,u) to 



a small but finite value in order to perform the integral in eq. (|6lJ) . This smoothens the numerically calculated free 
susceptibility (dot-dashed line) in comparison to the analytical expression given in appendix [b] (dotted line). 

Let us now study the imaginary part of the renormalized free susceptibility, Im X°(q , to), represented by the con- 
tinuous line in Fig. |] for half filling, U = 3 and q = ir/2. On the outer boundary ( |78^ it behaves essentially in the 
same way as Imx ilee (q, to), whereas the inner boundary ( f79| ) is completely washed out due to the renormalization of 
the occupation numbers rika- 



As will be discussed in section IV G , the strong coupling limit of our theory can be calculated analytically for half 
filling. Scaling our strong coupling result for lmx°(q,to) down to U — 3 leads to the dashed line in Fig. |] For U — 3, 
which is lower than the bandwidth and has thus to be considered as an intermediate interaction, Im^ of the complete 
renormalized RPA calculation (continuous line) qualitatively already resembles the properly scaled strong coupling 
result (dashed line). 

Fig. a shows Imx°(g,w) for U = 6 and half filling. Comparing the continuous with the dashed line, we notice that 
the renormalized RPA result for Im \° now agrees also quantitatively very well with the corresponding strong coupling 
result. The remaining difference comes from the finite imaginary part "i0 + " used for the numerical computation of 
Imx°(<7,w). 



2. Real part 

Let us now examine the real part of the free susceptibility x lrc °(<li w )- As can be seen from the analytic expressions 
given in appendix Q, Rex flee (q,to) diverges at the two boundaries of Imx flcc (g, to) which demark the particle-hole 
continuum. 

The divergence on the upper boundary, given by eq. (|78|) , makes the denominator of the pure RPA charge suscep- 
tibility (see eq. (pi))), 1 — Ux hcc {q, to), vanish at an energy above the continuum limit where an undamped plasmon 
is created. 

This is illustrated in Fig. ^ by the intersection of the horizontal line at +1/U with the real part of x frcc (<Z, w ) 
at u » 3.3. The dotted and dot-dashed lines represent the analytical or numerical expressions for Rex e (<7)<^)> 
respectively. Hence, as U is increased, the horizontal line at l/U is lowered and the pure RPA plasmon is shifted 
towards higher energies. 
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A similar scenario can be established for the longitudinal spin response in the pure RPA. Its denominator, 1 + 
17% vanishes at an energy below the lower continuum limit, given by eq. (|79|) , and an undamped magnon is 

created. Again, this is shown in Fig. |^ by the intersection of the horizontal line at —1/U with Rex frcc atwsi 1.6. 
As the interaction is increased, the magnon is shifted towards lower frequencies, and as it reaches u) = 0, the system 
becomes unstable. This will be discussed in more detail below. Note that slightly below the upper continuum edge, 
the dotted line representing Rex free also meets the horizontal line at —1/U. However, as Imx frcc is large in this 
region near the square-root singularity at the upper continuum boundary, it will not contribute to the pure RPA spin 
response. 

Fig. and Fig. || display the position of the plasmons and the magnons in the (q, w)-plane for U = 3 and for 
half and quarter filling, respectively. Once more, we notice that the pure RPA plasmons lie above the particle hole 
continuum, represented by the dotted area. The pure RPA magnons occur below the lower continuum boundary. As 
the lower boundary goes to u) — for q — > 2kp, the spin pole will meet the momentum axis near 2kp producing the 
well-known Peierls instability. This instability occurs in the pure RPA for infinitesimal U exactly at q = 2k p. When 
increasing the interaction, the region of instability is enlarged covering a growing interval around q = 2kp. 

It extends to the whole g-axis when U exceeds some critical interaction U stoneT , provided by Stoner's mean-field 
theoryEil. In mean field, the paramagnetic state, m — {n-\) — {ni) = 0, is always a local extremum in the energy 
surface as a function of the magnetization m. This extremum is a minimum for 

u < ^Stoncr = 27TSin^, (80) 

with n denoting the number of electrons per site. For U > U Stoacl it is a maximum. Nevertheless, it should be 
remembered that there is a range of interactions JJ MF < U < U SioneI for which the paramagnetic state is a local but 
not the global minimum, and that the global minimum is reached in the fully ferromagnetic state n = m. 

We thus conclude that the pure RPA only produces valid results in the region where no instabilities occur. Calcu- 
lating quantities having a contribution from the unstable region around q = 2kp, like e.g. the occupation numbers, 
does not make sense even for small interactions. Above U StoneT , the pure RPA is unstable for all q and ui. As we will 
see below, it is the virtue of the renormalized RPA (and also the SCRPA) to cure these instabilities, rendering possi ble 



the computation of the occupation numbers from eq. (|53|). These arguments will be underlined in section IV F by 
considering the energy weighted sum rule. 

Above the outer continuum limit (|78|), the real part of the renormalized free susceptibility, Rex°(g,w), behaves 
qualitatively like Rex (q, cu). In the same manner as above, we find an undamped plasmon at the energy where the 
denominator of the renormalized charge susceptibility in eq. (|64|), 1 — Ux°(q, w), vanishes. For small U, the plasmon 
is only slightly shifted to lower energies with respect to the pure RPA plasmon. For larger U, the renormalization 



effects arc stronger. We will see in section [V G that the frequencies of the renormalized plasmons remain finite as U 
goes to infinity. The pure RPA plasmons, in contrast, occur at infinite frequencies in the limit of large U. 

Below the outer continuum limit (|78|), the renormalization effects in Rex°(9, oj) are more drastic. This is shown by 
the continuous line in Fig. |^ for U — 3 and half filling. The singularity of Re ^ frcc at the lower continuum boundary 
(u ss 2) is completely damped and Rex° is an almost structure-less function within the renormalized continuum, i.e. 
from uj = up to almost the continuum edge (|7§| ) at uj ss 2.8. 

This gives, even for small interactions, rise to qualitative changes in the longitudinal spin response. At the energy 
where the denominator of the longitudinal spin susceptibility, 1 + J7x°(g, w), becomes resonant, the imaginary part 
of x°(?) w ) is large. Therefore, the renormalized spin response will show a broad maximum instead of the undamped 
magnon found in the pure RPA spin response. 

Apart from the qualitative differences expected between the pure and the renormalized RPA spin responses, there 
are important consequences for the regime where the pure RPA is unstable. We recall that in this g-range, the pure 
RPA magnon frequency becomes imagi nary and the magnon pole in the pure RPA spin response disappears. As one 



consequence, we will outline in section IV F that the energy weighted sum rule for the pure RPA spin response will 
be violated. 

In contrast, the broad maximum in the renormalized RPA spin response will persist even for the q— vectors where 
the pure RPA is unstable. Therefore, the renormalized RPA spin response fulfills the ener gy w eighted sum rule for 



every q, even for interactions U > fj stoncl ' . This will be explained in more detail in section IV F. By this means, the 
Peierls instability occurring in the pure RPA is cured by the renormalization process. 

In the limit of large interactions, the real part of x° becomes completely flat within the renormalized particle- 
hole continuum, and 1 + £/Rex°(<7, uj) vanishes within the whole continuum. As Imx°(q,oj) is smallest within the 
continuum for uj — > 0, the broad maxi mum in the longitudinal spin response is shifted towards zero frequency. The 



consequences are discussed in chapter [V E 



We finally remark that even for small U, the renormalization effects arc strong enough to lock Rex ^,^) to values 
greater or equal —1/U (see Fig. ||). Consequently, the scenario given above is valid, and the renormalized spin response 
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shows already for weak interaction strengths a broad continuum peak rather than a sharp magnon pole. 



D. Charge response 

Due to the rather small renormalization effects of u>) outside the particle-hole continuum, we expect the charge 
response in the pure and in the renormalized RPA to behave similarly. 

For intermediate interactions, the pure RPA charge response is given by a rather small continuum which is limited by 
the two boundary conditions ([78]) and (|7^) , and an undamped plasmon lying above the continuum. This is illustrated 
for U — 3 and half filling by the dotted line in Fig. ||, where we used the analytical expressions for x frc °j given in 
appendix |b| to evaluate the pure RPA charge res ponse . Evaluating the pure RPA charge response numerically, i.e. 



with a finite value of "i0 + " as discussed in section [V A , leads to the dot-dashed line in Fig. g. We observe that, due 
to the finite "i0 + ", the pure RPA plasmon at u) « 3.3 is broadened with respect to the analytical curve, and that the 
sharp cutoff at the continuum boundaries is smoothened. 

Comparing the continuous line with the dot-dashed line in Fig. ^| sho ws that the renormalized charge response 



agrees essentially with the pure RPA predictions. As explained in section IV C 2, the main differences are that the 
tail of the renormalized charge continuum now goes down to u> — and that the plasmon is slightly shifted towards 
lower energies. The precise position of the renormalized plasmon is indicated in Fig. |^ by a continuous vertical line 
in the center of the numerically computed plasmon pole which, for numerical reasons, has a finite width. 

Qualitatively, the charge response for U = 3 is already very similar to the properly scaled strong coupling limit of 
our theory. This can be seen by comparing the dashed and the continuous lines in Fig. Moreover, the plasmon 
position of the pure RPA and the strong coupling limit of the renormalized RPA agree very well, such that the dotted 
and the dashed vertical lines cannot be resolved from another. Note, however, that this agreement is purely accidental. 

The renormalization effects become more drastic as the interaction strength is increased. For half filling and U = 6 
Fig. |l0| illustrates that the renormalized charge continuum now not only goes down to u> = (continuous line) but also 
is much stronger than the pure RPA charge continuum (dot-dashed line). Further, the energy shift of the plasmon 
obtained from renormalization is larger than in the weak coupling limit. 

For U = 6, the charge response already agrees quantitatively very well with the strong coupling limit of our theory. 
Comparing the continuous and the dashed line sh ows th at both, the continuum contributions and the positions of the 



plasmon peaks are in good agreement. In section IV G we will discuss in detail that the renormalization prevents the 
plasmon pole in the large U limit from being shifted to infinite frequencies as this happens for the pure RPA plasmon. 
Similar results are obtained away from half filling. 



E. Longitudinal spin response 



According to the discussion of the real parts of x free (<?> w ) an d X°(<7> w ) m section 1VC 2 , we expect the spin response 
to change qualitatively when passing from pure to renormalized RPA, even in the regime of weak and intermediate 
coupling. 

This is illustrated in Fig. [ll] for q = tt/2, {7 = 3 and half filling. For these values there are already significant 
changes, although the pure RPA is still stable. The pure RPA spin response (dot-dashed line) consists of a rather 
small continuum being restricted to the area in between the two boundaries given by eqs. ( |7q ) and (|79|), and an 
undamped magnon occurring below this continuum. 

As explained before, the broadening of the magnon is of numerical origin. The corresponding analytic expression 



can be obtained using the representation of x flcc from appendix 0. The latter leads to the dotted line in Fig. 1 1 



It was argued in section IV C 2 that the magnon pole disappears during the renormalization procedure. Hence, 
the renormalized RPA spin response is fully described by a continuum exhibiting a broad peak (continuous line). 
Like the renormalized charge continuum, this spin continuum starts at u> = and goes up to the upper continuum 
boundary ([78]). The dashed line in Fig. [ll] shows the spin response in the strong coupling limit of our theory scaled 
to U = 3. A comparison with the renormalized RPA result (continuous line) shows that especially for low frequencies 
there are still important differences. Thus, at U = 3 the strong coupling limit is not yet reached. 

If U is increased or if the g-vector is chosen in the domain around 2kp where the pure RPA is unstable, the pure 
RPA magnon frequency will become purely imaginary. Fig. [l^ shows this case for U = 6, q = tt/2 and half filling. The 
pure RPA spin response is represented by the dot-dashed and the underlying dotted line, depending on whether the 
numerical or the analytical expression is monitored. We notice that the magnon peak in the pure RPA spin response 



vanishes completely. This corresponds to an unphysical situation, as will be underlined in section IV F by sum rule 
arguments. 



18 



In renormalized RPA, the broad maximum in the longitudinal spin response is shifted towards zero frequency 
when the interaction is increased. This can be seen by comparing the continuous lines in Fig. a nd Fig. [l2], which 
represent the renormalized spin response for U — 3 and U = 6, respectively. Whereas in FigTo, the renormalized 
spin response reaches its maximum at a finite frequency we see from Fig. |l2j that for U = 6 it is strongly peaked at 
lu = 0. This behaviour is characteristic for the strong coupling limit of our theory, given by the dashed line which 
cannot be resolved from the continuous line in the latter graph. In renormalized RPA, the energy weighted su m rule 
is fulfilled for all U, as will be demonstrated in the next subsection. Nevertheless, we will point out in section [V G 
that the singularity of the spin response at to = may lead to divergences in correlation functions. 



F. Energy weighted sum rule 



In section [II B a sum rule is derived for the energy weighted spin and charge response. In the one-dimensional 
Hubbard model, eq. ( |68| ) connects an energy weighted integral over the response functions to the mean kinetic energy 
per site, ( i ), times a form factor: 

00 

Sf{q): -2 (t) (1-cosg) = -- J duwlm^fcw) 

— OO 
OO 

S?(q): ~\ (i) (1-cosq) = ~ f d^Im^( ? , U ). (81) 



In section [II B it is argued that it not only holds true for the exact response functions x ch (<Z, w) and x sp (<7, w) but 
as well for the response functions in Self-Consistent and in renormalized RPA if the mean kinetic energy per site, 
(i), is calculated self-consistently from the corresponding Green's function. 

Moreover, it is well known that the sum rule for the pure RPA resp onse functions is fulfilled if ( t ) is calculated 
with the Hartree-Fock ground-state as long as no instability occMral3'E3. 

In Fig. |lj, [l5| and |l^ we show the result of the sum rule checks for U = 3, and for half and quarter filling, 
respectively. In all figures, the left hand sides of the eqs. ( |3l| ) are plotted with the dotted lines. As their (1 — cosq)- 
behaviour is independent of the approximations made in the Green's function, we will consider them as "reference 
lines". Nevertheless, they are scaled with a prefactor, (i), which depends on the Green's function. In renormalized 
RPA, the mean kinetic energy per site is less negative than in Hartree-Fock, since the momentum distribution is 
smoothened. Thus, the reference lines for the renormalized RPA sum rules lie always below the pure RPA lines. 

For the renormalized RPA susceptibilities, the right hand sides of the eqs. (|8l]) are represented by the continuous 
lines. They cannot be resolved from the corresponding dotted reference lines. Hence, in renormalized RPA the sum 
rule is fulfilled for all g-vectors and for both, the charge and the longitudinal spin response. 

Calculating the rhs. of the eqs. ( |8l| ) with the pure RPA susceptibilities yields the dashed lines. The sum rule for the 
charge response is monitored in the figures |l3| and [l5|. There, the dashed lines cannot be resolved from their dotted 
reference lines. This means, that the sum rule for the pure RPA charge response holds true for all g-vectors. 

In the figures [l4| and [ll| we show the sum rule check for the pure RPA spin response. The dashed lines, which 
represent the energy weighted pure RPA spin response, agree with their reference lines, apart from a range of g-vectors 
around the Peierls vector 2k p. This q-range coincides with the g-range where the corresponding magnon dispersion 
in Fig. |t| and Fig. |[ respectively, goes to ui = 0. Therefore, the energy weighted sum rule for the pure RPA spin 
susceptibility is only fulfilled for the range of momenta where the pure RPA is stable. 

Again, we want to emphasize the quality of the renormalized and the Self-Consistent RPA to restore the /-sum 
rule for all momenta and interactions. 

G. The large U limit 

At half filling and for large U , the occupation numbers predicted by renormalized RPA can be fitted very accurately 



by 



n ka = i M + ^cosfc j (S2) 
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whereas the large U expansion of the Bethe ansatz solution yield; 

1 ( 81n2 \ 
n k a = g I 1 H — cosfc J . (83) 

Comparing the dotted line and the lowest continuous line in Fig. || shows that for U — 5, the expression (|8^) agrees 
already very well with the numerical results for nu- We will show in the following that eq. (^2|) indeed is the fully 
self-consistent solution of the renormalized RPA equations in the strong coupling limit. 

Based on eq. (p^), it is possible to give an explicit expression for the renormalized free susceptibility. As the con- 
tinuum boundary condition for step-like occupation numbers, eq. (f79|), becomes meaningless in the limit of strong 
interactions, the only characteristic energy is provided by the upper continuum limit ([78]). Hence, the explicit expres- 
sions for X (<?, w) do not depend on q and u> as independent variables anymore, but can be denoted as a function of 
one single variable £, which is the energy in units of the upper continuum limit (JTS 



|4sin f | 



(84) 



Performing the integration in eq. ( |6l| ) with occupation numbers as given by eq. ( |S2| ) yields the renormalized free 
susceptibility 

Like for the free susceptibility x frcc , the imaginary part of x° still has a square root singularity at the continuum limit 
|£| = 1, whereas the real part of x° is completely flat within the particle-hole continuum. By scaling the renormalized 
free susceptibility of the strong coupling limit, eq. (|8q), to U = 3 or U = 6 we obtain the dashed lines in Figs. |[ || 
and ||. 

Substituting eq. j8^) in eq. ([m] ) yields the renormalized RPA response functions. The large U charge response is 
characterized by a small continuum within |£| < 1, 



ImX ch (e) = -| • (86) 

Outside the particle-hole continuum, the charge response has a collective pole at £ = 2/ v3, describing a undamped 
plasmon obeying the dispersion relation 



Qp V3 



. q 

sin — 

2 



(87) 



Notice that the dispersion relation of the renormalized plasmon becomes independent of U in the large U limit. The 
pure RPA, in contrast, predicts in this limit a plasmon dispersion proportional to U . The strong coupling charge 
response scaled to the appropriate U is represented in Fig. ^ and Fig. |l^ by the dashed lines. 

From the second of the eqs. (^34|), in combination with the renormalized free susceptibility J85|), we obtain the large 
U spin response in renormalized RPA. It is fully described by a continuum for |£| < 1, 



Im X sp (0 = ~ , (88) 

and vanishes elsewhere. This is plotted with dashed lines in Fig. |ll| for U = 3 and Fig. |lj for U — 6. No collective 
excitations occur, since the denominator of the renormalized RPA spin response, l+U~x°, is always finite. Nevertheless, 
we see from eq. ( p8[ ) that the spin continuum diverges at to — with consequences that we will discuss below. 

It can be shown that the real and imaginary part of x° from eq. ( |S5|) fulfill the Kramers-Kronig relations. Moreover, 
we will show in appendix ^ that the occupation numbers calculated from the large U expression of the renormalized 
RPA Green's function via the spectral theorem (E^) are consistent with eq. j82|). This means that the occupation 
numbers (|83) fulfill together with the strong coupling susceptibilities fl86| ) and (B8[) the self-consistency condition of 
the renormalized RPA. 

The number of double occupancies per site are given by 
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jf^n^nn = (n T ) + — , (89) 

i 

where u> c stands for the correlated potential energy per site, 

WC = JpYsi a t] a k+q1 a t+qi a Pl ) ■ ( 90 ) 
kpq 

w c is independent of U, since the large U Green's function scales with 1/U. 

Nevertheless, as the expectation value in eq. (Q) can be obtained from the spectral theorem fl54| ) by integrating 
the renormalized RPA Green's function over positive frequencies, the 1/ui divergence of the spin response (|8|) makes 
the expectation value w c diverge logarithmically. This means that both, the total potential energy per site and the 
number of doubly occupied sites are going to — oo as U — > oo. 

One way to correct these deficient results is-,to calculate the ground-state energy E(U) from the kinetic energy 



T(U), using the Hellmann-Feynman-Theoremc 



oo 

E(U) = uJdy^. (91) 

u 

The behaviour of the kinetic energy in the renormalized RPA can be obtained from the occupation numbers js^). 
For large U, this yields 

T(U) = -1 . (92) 



Again, only the prefactor differs from the exact result, 



•(to - ~ , m 



known from the large U expansion of the Bethe ansatz solutionsl 
Calculating the ground-state by eq. ([n]), our theory predicts 

E(U) = -i (94) 

as U — > oo. The large U expansion of the renormalized RPA ground-state energy of the half filled Hubbard chain is 
shown in Fig. [l^by the dashed line. For arbitrary interactions, the renormalized RPA ground-state energy calculated 
with the Hcllmann Fcynman theorem is illustrated by the "x" -symbols. This compares reasonably well with the 
exact result, known from Bethe ansatz (continuous line) and its large U expansion (dot-dashed line). 

We now have access to the potential energy as the difference between ground-state energy and kinetic energy. It 
thus behaves as 2/U for large U. This implies that the number of double occupancies vanishes as 2/U 2 for strong 
interactions. In Fig. [l8|, the number of double occupancies at half filling is monitored as a function of U. We see that 
for small interactions, the renormalized RPA ("x" -symbols) reproduces the exact result (continuous line). For large 
interaction strengths, our theory matches the Bethe ansatz expansion, — (41n2)/C/ 2 (dot-dashed line), apart from the 
prefactor. 



V. DISCUSSION, CONCLUSIONS AND OUTLOOK 

In this work, we performed a first application of the Self-Consistent RPA (SCRPA) theory to the Hubbard model. 
In itself, the SCRPA is an approximation to the general Dyson Equation Approach (DEA) to correlation functions 
where the full (exact) mass operator is replaced by its instantaneous contribution. For the single particle Green's 
function this strategy leads to the standard Hartree-Fock theory. In analogy, it has been argued in the past that the 
SCRPA corresponds to a HF theory for fermiorupair clusters. Recently, this theory has produced very interesting 
results in various domains of many-body physics&Q. 

Unfortunately, being a (non-linear) mean-field theory for non-local correlation functions, SCRPA is numerically 
very demanding. As a first step, we therefore had to proceed to further rather drastic simplifications. The latter 
consist in retaining correlations only in the single-particle occupation numbers. This approximation to SCRPA is 



21 



known in the literature as renormalized RPAiB In spite of this, the essentials of the self-consistency and closure 
aspects remain intact. As a further virtue, the /-sum rule is shown to be fulfilled in Self-Consistent as well as in 
renormalized RPA. This also implies that the Goldstone theorem is fulfilled and symmetries are treated correctly. 

We solved numerically the renormalized RPA equations for the one-dimensional single-band Hubbard model in 
the paramagnetic phase, for different fillings and interactions. Although we were aware of the difficulty of describing 
one dimensional models because of the extreme importance of quantum correlationsu, there were multiple reasons 
for this choice. In a first place, the exact solution of Lieb and W ueI provides a benchmark for our results, which, in 
higher dimensions, does not exist. As was argued in section [II A, we expect the renormalized RPA to perform better 
with increasing dimensionality. Therefore, Id can be considered as a "worst case check" for our approximation. The 
second reason is mainly technical: The self-consistency equation for the occupation numbers (^3|) illustrates that the 
numerical effort increases with the square of the spatial dimension. Therefore, the experience in Id is desirable before 
attacking higher dimensions. In the last place, even in one dimension, we are able to test explicitly the essentials of 
our method and the convergence of the iteration cycle. 

As expected, most of the particularities of the one-dimensional model, as e.g. Luttinger liquid behaviour away from 
half filling, could not be reproduced. The scenario predicted by our theory rather confirms the results wjhich, were 
obtained from methods designed for higher dimensional models, like e.g. the Hubbard-III approximatiorJIS'Lj: For 
small U, we find a Fermi-liquid like metallic ground state. The system undergoes a Mott-Hubbard transition for an 
interaction slightly smaller than the bandwidth. The exact value of the critical interaction depends on the filling. For 
stronger interactions, our theory predicts an insulator for all fillings. 

Despite these deficiencies, we obtain several interesting results. In the strong coupling regime of the half filled 
model, for example, we are able to express the central quantities of our theory analytically. The renormalized RPA 
momentum distribution function, given by eq. (|82]), agrees, apart from a prefactor, with the oc cosfc behaviour 
known from the large U expansion of the Bethe ansatz solution. 

Moreover, the mean-field spin instability around 2fcp, which causes the breakdown of the pure RPA for any finite 
interaction, turns out to be cured in renormalized RPA in the sense that no more purely imaginary eigenfrequencies 
occur in the spin channel. On the other hand, the renormalization is still rather weak such that a strongly overdamped 
spin pole remains at low energy. These low-lying spin excitations give rise to a slow (logarithmic) divergence of the 
two-body correlation functions, which is carried into the number of double occupancies and thus also into the ground 
state energy. The stronger renormalization contained in the SCRPA would certainly cure this pathology, since it also 
renormalizes (screens) the interaction self-consistently. In this sense, the SCRPA bears some similarities with the 
approach of Tremblay et al.u. 

Several of the renormalized RPA results may nonetheless be improved by applying the Hellmann-Feynman theorem. 
Using this approach, which may simulate a step towards the full solution of the SCRPA, the ground-state energy 
compares for all interaction strengths reasonably well with the exact results, known from Bethe ansatz. Moreover, 
the number of doubly occupied sites shows a qualitatively correct behaviour over the whole U range. Especially for 
U — ► 0, the mean-field value is recovered, and, for U — > oo, the double occupancies vanish like l/U 2 , as predicted by 
the exact solution. 

In the strong coupling regime of the half filled model, this approach reproduces the exact results for the momentum 
distribution function, the ground-state energy, the kinetic and potential energy, and the number of doubly occupied 
sites, apart from a general prefactor. We would obtain the right prefactor by substituting the bare Hubbard U with 
a screened interaction, or, in other words, by multiplying U in our theory by a factor 1/(2 In 2) w 0.72. 

Since in renormalized RPA the Hubbard interaction is approximated by a simple spin-flip interaction, it is not 
astonishing that the best results are produced at half fillingcj. Away from half filling, other scattering processes, still 
taken into account in SCRPA but neglected in renormalized RPA, become important. Therefore, the renormalized 
REA is less effective. The Luttinger liquid behaviour, exhibited by the exact momentum distribution even for U — > 
ooc3, is not obtained. In the strong coupling limit, the renormalized RPA produces smooth occupation numbers, 
although still having a steeper slope at kp as at half filling. For small interactions, the renormalized RPA still 
predicts a discontinuity in the momentum distribution at kp. This implies that the Luttinger theorem is satisfied. 

A further point to be discussed is the fact that we completely neglected the dynamical interaction part, %a%(v\> 
defined in eq. (11). It is certainly true, as recently pointed out by Logan for the infinite dimensional Hubbard modelES, 
that this dynamical interaction should be taken into account, in order to describe the low energy scale in the spin 
channel correctly. Nevertheless, in the present work, some dynamical effects in the spin channel are considered via the 
coupling of the occupation numbers to the RPA ground state correlations. A full inclusion of the dynamical effects 
goes, however, beyond the scope of this paper. 

At this point, it may be appropriate to shortly come back to some technical aspects of the SCRPA which we did 
not develop in the main text in order not to overduely extend the size of the paper. In the introduction, we mentioned 
that the Equation of Motion Method (EOM), on which our formalism is based, goes back rather far in time. The 
first major theoretical input was developed by D.J. Rowe in his review articlea, where the calculation of density- 
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density correlation functions-is described in the context of nuclear physics. The method was later applied to strongly 
correlated electrons by Rotbu. She evaluated the single-particle Gceen's function by coupling it in an approximate 
way to the three particle propagator (see also Beenen and EdwardstZl for a more recent application to the Hubbard 
model). Since then, in solid state physics, the EOM has, to the best of our knowledge, exclusively been used for the 
calculation of single particle properties!^. 

As we point out in ref.u, the optimal procedure will be to combine single-particle and fermion-pair channels. 
Indeed, as is well known, the single particle mass operator can be expressed exactly by the two-particle T-matrixKJ. 
Replacing it by its SCRPA counterpart then naturally leads to a self-consistent coupling of both channels. In this 
scheme, as we will describe in more detail in a future publication, the single-particle occupation numbers will not be 
evaluated from the particle-hole propagator (see eqs. (g7j) and (|53|), respectively) but directly from the single-particle 
Green's function. A further advantage of this strategy is that single-particle and fermion-pair properties are obtained 
simultaneously and on equal footing. We did not follow this route in this paper, since it again would have strongly 
increased the numerical difficulties. In spite of these possible improvements, the present investigation shows that 
the essentials of the SCRPA theory (i.e. the self consistency procedure) work correctly in a numerical application 
to a homogeneous system of strongly interacting fermions. Once the method will be solvable in its full complexity, 
the possible applications are very numerous. Indeed, the SCRPA is a very flexible formalism applicable to strongly 
correlated Fermi systems, but also to Bose or spin systems. The study of such systems are planned in the future. 

A very appealing application of the SCRPA may be the Hubbard model in infinite dimensions, since, on one hand, 
the physics in d = oo is expected to be somewhat similar to d = 3. On the other hand, spatial fluctuations are 
suppressed in d = oo which allows to reduce the many-body problem to either a dynamical single-site problem, or, 
to an effective one dimensional problem. We thus may expect, that the numerical solution of the SCRPA will be 
feasible in infinite dimensions, in contrast to finite dimensions, where the effort will grow with an exponent 3d, since 
the effective Hamiltonian contains correlation function depending on three momenta. 
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APPENDIX A: DERIVATION OF THE SCRPA FROM A VARIATIONAL PRINCIPLE 



In this section, we will briefly outline the derivation pof the particle-hole SCRPA equations from a variational 
principle. An analogous method was derived by Barangeiuj for the single particle case. 
Let us therefore consider the spectral representation of the retarded Green's function 



E 



X\n)(n\X+\0) (Q\X+\n){n\X\0) 



UJ - UJ n Q + i0+ 



UJ + UJ n0 + i0 + 



(Al) 



where |0) is the exact ground state of H, and uj n o denotes the excitation energy for the exact eigenstates |n). For the 
particle- hole problem, we set e = +1. The excitation operators X + are given by 



x + = 22 x kp a i a p ■ 
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The normalized mean excitation energy is given by: 

J> n0 (\(n\X+ 

Si = 



- \(n\X 



E(\(n\X+\0)\ 2 - \(n\X\0)\' 

n ^ 

oo 

J dojujIm((X;X+))l ct 



J dujlm((X;X+)) r u 



(A3) 



23 



The equivalence between the first and the second line can be seen by substituting the spectral representation (Al) 
in the second line. The denominator can easily be evaluated with the spectral theorem ([lB]). This yields the norm 
(0| [X, X + ] |0), which wo uld b e equal to unity if X + were ideal bose operators. 

We now minimize eq. (A3) with respect to the excitation operators. The X + with the lowest mean excitation 
energy obey the condition: 



dSi 
dx kp 







(A4) 



It is straightforward to verify that eq. ( A4) leads to the SCRPA for the particle- hole propagator derived in section [I C 



APPENDIX B: FREE SUSCEPTIBILITY IN ONE DIMENSION 



The free susceptibility is obtained by performing the integral in eq. (fn]) for step-like occupation numbers J63j). 
Because of the symmetry of the dispersion relation the real part of X\Qi u ) is symmetric in q and u>, whereas the 
imaginary part is symmetric in q and antisymmetric in to. 

In the paramagnetic phase, = n^i, and in one dimension, the explicit expression for the real part is: 
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with z = (4 sin |) 2 — J 1 and 



arctan x for z < 
atanx = ^ artanhx for z > and |x| < 1 
I arcothx for z > and \x\ > 1 



For the imaginary part, we find in agreement with Benard et al.0 
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APPENDIX C: SELF-CONSISTENCY OF THE LARGE U LIMIT 

In this section, we will briefly outline that the occupation numbers which are given by eq. (|8^) for the half filled 
Hubbard chain in the large U limit are indeed a fully self-consistent solution of the renormalized RPA equations. 

If we assume the occupation numbers from eq. J82|), we find the renormalized free susceptibility (^5|) by calculating 
the fc-sum in eq^J6l|). The large U limit of the renormalized RPA Green's functions is then obtained by substituting 



X°(?,w) in eq. (|2|) 

We are now able to calculate a new set of occupation numbers by inserting this renormalized RPA Green's function 
in eq. (|53]). For convenience, however, we will rather use the commutator spectral theorem, eq. (|l5|), itself: 
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J Im g kcrpa (q,uj) 

(Uktr - n k- 



div Im 



xold 



w - [£ fe+g -£fc] + ^0+ l-?7x?(?,w)?7x;(g,w) 



(CI) 



For the self-consistency to be fulfilled, we now have to show that the "new" occupation numbers are equal to the 
"old" ones. As we are in the paramagnetic phase, we may drop the spin indices and convert the integrand in a partial 



fraction. After substituting x 
occupation numbers: 



;/ 4sin | and y kq = [sk+ q — £&] / 4sin | , we find for the ratio between new and old 



oo 

I(Vkq) = --^ J dxlm 



x - y kg + i0+ \1-U X ° 1 + U X ° 



(C2) 



In the following, we will treat the two fractions in the integrand separately, considering in analogy to eq. ( p4| ) the 
first term as the charge term, and the second term as the spin term. With the renormalized free susceptibility from 
eq. (p5D, the denominator of the charge term may be written as 



1 - U X °(x) 




In the same way, the denominator of the spin term is 

-i 



l + U X Q (x) = 



for |a;| < 1 
+ «0 + signa; for \x\ > 1 

for Id < 1 



, 13 — iO + signa; for \x\ > 1 

yx 2 — l 11 



(C3) 



(C4) 



To solve the integral, we have to account for three different contributions: The first one comes from the pole of 
the free ph Green's function, l/(x — ykq + «0 + ), which lies always in the ph continuum, i.e. \ykq\ < 1- With Dirac's 
identity, 



1 



we find for the charge contribution 



f{x) ± iOH 



ifiVkq) 



= V 



1 



fix) 



T ivS [f(x)] 



Via 



tyiq 



(C5) 



(C6) 



The corresponding spin contribution is zero, since 1 + Ux°(x) is purely imaginary within the ph continuum. 

Secondly, we have to integrate over the charge and spin continuum, respectively. Therefore, we combine the 
imaginary part of 1/(1 ± Ux°) with the real part of the free ph Green's function. This yields 



itiVkq) 



3<4 



6 4 



3*2.) 



(C7) 



The last contribution comes from the collective poles. Again, the spin contribution is zero, as th e reno rmalized spin 
response does not show any magnon peak in the stro ng c oupling limit. As pointed out in section IV G , the plasmon 
peak occurs at x = ±2/v3< Using the Dirac identity JC5D once more, we get 



rch 



(Vkq) 



3 4 - 3yi 



(C8) 



Summing the three contributions, (|C6|), ( |C7|) and ( |C8| ), yields unity for the ratio of the new and the old occupation 
numbers, I(yk q ). By this means, we have found a fully self-consistent solution of the renormalized RPA equations. 
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FIG. 4. Imaginary part of the free and renormalized free susceptibility x (<7 = n /^i u} ) f° r the half filled Hubbard chain at 
[7 = 3. 
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FIG. 6. Real part of the free and renormalized free susceptibility x°{q = t/2, w) for the half filled Hubbard chain at U — 3. 
The intersection with the horizontal lines at ±1/U indicate where the collective excitations occur (see text). 
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FIG. 7. Plasmon and magnon dispersion for the half filled Hubbard chain at U = 3. The renormalized RPA plasmons are 
plotted with circles. "+"- and "x"-symbols stand for the pure RPA plasmons and magnons, respectively. The dotted area 
illustrates the non-interacting particle-hole continuum. Its boundaries are given by the continuous lines. 
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FIG. 9. Charge response x ch (<l ~ nfti w ) m pure and renormalized RPA for the half filled Hubbard chain at U = 3. 
The vertical lines indicate the plasmon peaks. Accidentally, the pure RPA plasmon and the strong coupling plasmon occur at 
almost the same energy, such that the dotted and the dashed vertical line cannot be resolved from another. The thin continuous 
vertical line illustrates the precise position of the renormalized RPA plasmon, which itself is represented by the singularity in 
the renormalized RPA charge response (continuous line). 
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FIG. 10. Charge response x ch (<3 = n/2,u>) in pure and renormalized RPA for the half filled Hubbard chain at U = 6. The 
vertical lines indicate the plasmon peaks (see Fig. . The strong coupling limit of our theory is reached, and the renormalized 
RPA plasmon (continuous vertical line) cannot be resolved from the strong coupling plasmon (dashed vertical line). 
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FIG. 13. Energy weighted sum rule for the renormalized and pure RPA charge susceptibility for the half filled Hubbard 
chain at U = 3. The dashed and continuous line correspond to the rhs. of the charge sum rule Si h (q), eq. ([H]), computed 
with the pure and renormalized RPA charge susceptibility, respectively. The left hand side of the sum rule is plotted with 
the thin dotted "reference" lines. As the sum rule is fulfilled in both cases, these reference lines cannot be resolved from the 
corresponding continuous or dashed lines. 
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FIG. 14. Energy weighted sum rule for the longitudinal spin susceptibility in renormalized and in pure RPA for the half 
filled Hubbard chain at U = 3. The dashed and continuous line correspond to the rhs. of the spin sum rule Sl P (q), eq. (jil|), 
computed with the pure and renormalized RPA spin susceptibility, respectively. The left hand side of the sum rule is plotted 
with the thin dotted "reference" lines. As the sum rule for the renormalized RPA is fulfilled, the reference line cannot be 
resolved from the continuous line. For the pure RPA, there is a region around 2&f where the sum rule breaks down. 
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FIG. 15. Energy weighted sum rule for the renormalized and pure RPA charge susceptibility for the quarter filled Hubbard 
chain at U = 3. Description as in Fig. Il3. 



ren. RPA, spin 




FIG. 16. Energy weighted sum rule for the longitudinal spin susceptibility in renormalized and in pure RPA for the quarter 
filled Hubbard chain at U = 3. Description as in Fig. [14l 
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FIG. 17. Ground state energy per site for the half filled Hubbard chain. Bethe ansatz expansion from Baeriswyl et al.@ 
Large U limit of the renormalized RPA from eq. (|94|). 
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FIG. 18. Number of double occupancies per site for the half filled Hubbard chain. 
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